Load required packages for entire analysis
dlnm.pkg <- c("tidyverse", "lubridate","xts","ggthemes","PerformanceAnalytics","reshape2","rugarch","timetk","parallel","timeSeries","tseries","data.table","ggplot2","dlnm","broom","caret","gridExtra","splines","splines2","pspline","cowplot","mgcv","spi","chron","gridGraphics","grid","pscl","MASS", "AER", "Hmisc", "MuMIn", "VGAM", "forecast", "seasonal", "plotly", "ggmap", "rgeos", "tmap", "maptools", "maps", "ggfortify", "htmltools","webshot","knitr","flexdashboard", "imager", "httr", "curl", "here")

lapply(dlnm.pkg, library, character.only=TRUE)
Load Malawi map dataset
dlnmtmp <- tempfile()
download.file("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/malawi_map.zip", destfile=dlnmtmp)
unzip(dlnmtmp, exdir = ".")
malawi.map <- rgdal::readOGR(".","malawi_map")
Subset Blantyre map
blantyre1.map <- malawi.map@data$OBJECTID >289 & malawi.map@data$OBJECTID <297 
blantyre2.map <- malawi.map@data$OBJECTID >308 & malawi.map@data$OBJECTID <311
blantyre3.map <- malawi.map@data$OBJECTID >342
Convert shape file to dataframe
blantyre.map <- rbind(fortify(malawi.map[blantyre1.map,]), fortify(malawi.map[blantyre2.map,]), fortify(malawi.map[blantyre3.map,]))
blantyre.map$id <- as.integer(blantyre.map$id)
Load demographics and features dataset
blantyre.demog <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_demog.csv"))
map.features <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_features.csv"))
blantyre.demog$id <- as.integer(blantyre.demog$id)
map.all <- merge(x=blantyre.map, y=blantyre.demog, by="id", x.all=TRUE)
rm(list = ls()[grep("^blantyre", ls())])
Plot the map of Blantyre
ggplot() + 
  geom_polygon(data=map.all,aes(x=long,y=lat,group=group,fill=popc),colour="gray50") + 
  theme_classic() + 
  theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) + 
  labs(fill="(1998 - 2008) Population censuses") + 
  xlab("Longitude") + 
  ylab("Latitude") + 
  geom_point(data=map.features, aes(x=long, y=lat, shape=Geolocation, size=Geolocation), color="black") +
  scale_shape_manual(values=c(17, 16, 3)) +
  scale_size_manual(values=c(2,4,3)) + 
  theme(legend.key.height=unit(0.8, "line")) + 
  theme(legend.key.width=unit(0.8, "line"))

Load case dataset and subset by serovar
case <-read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/case.csv"))
case$case_date <- dmy(case$case_date)
case$case_count <- c(1)
case.typhi <- subset(case, organism=="typhi")
case.iNTS <- subset(case, organism=="iNTS")
Assign 0 when no case was diagnosed
case.typhi <-aggregate(case.typhi$case_count, by=list(case.typhi$case_date), FUN=sum, na.rm=TRUE)
names(case.typhi) <- c("date", "case_count")
case.typhi <- merge(case.typhi, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.typhi[is.na(case.typhi)]<-0

case.iNTS <-aggregate(case.iNTS$case_count, by=list(case.iNTS$case_date), FUN=sum, na.rm=TRUE)
names(case.iNTS) <- c("date", "case_count")
case.iNTS <- merge(case.iNTS, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.iNTS[is.na(case.iNTS)]<-0
Convert to XTS and aggregate monthly
case.typhi = as.xts(case.typhi[,-1,drop=FALSE], order.by=as.Date(case.typhi[,1]))
case.iNTS = as.xts(case.iNTS[,-1,drop=FALSE], order.by=as.Date(case.iNTS[,1]))
case.typhi <- apply.monthly(case.typhi, FUN=sum)
case.iNTS <- apply.monthly(case.iNTS, FUN=sum)
Load climate dataset
climate <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/climate.csv"))
Define rain/temp between chileka & chichiri stations
climate$climate_date <- dmy(climate$date)
climate$rainfall <- (climate$chil_r+climate$chic_r)/2 
climate$temperature <- ((climate$chil_mint+climate$chil_maxt)/2+(climate$chic_mint+climate$chic_maxt)/2)/2
climate <- subset(climate, select=c(climate_date, rainfall, temperature))
climate.rain <- subset(climate, select=c(climate_date, rainfall))
climate.temp <- subset(climate, select=c(climate_date, temperature))
Convert to XTS and aggregate monthly
climate.rain <- as.xts(climate.rain[,-1, drop=FALSE], order.by=as.Date(climate.rain[,1]))
climate.temp <- as.xts(climate.temp[,-1, drop=FALSE], order.by=as.Date(climate.temp[,1]))
climate.rain <- apply.monthly(climate.rain, FUN=mean)
climate.temp <- apply.monthly(climate.temp, FUN=mean)
Create tibbles for time series decomposition
case.typhi <- tk_tbl(case.typhi, preserve_index=TRUE, rename_index="date") 
case.iNTS <- tk_tbl(case.iNTS, preserve_index=TRUE, rename_index="date") 
climate.rain <- tk_tbl(climate.rain, preserve_index=TRUE, rename_index="date") 
climate.temp <- tk_tbl(climate.temp, preserve_index=TRUE, rename_index="date")
Seasonal-adjusted NTS cases
case.iNTS.ts <- ts(na.omit(case.iNTS$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.iNTS.ts)))), preserve_index=FALSE)
case.iNTS$case_count_sea <- trend_n$value 
case.iNTS$case_count_sea[case.iNTS$case_count_sea < 0] <- 0
setnames(case.iNTS, old="case_count", new="case_count_obs")
trend_n <-tk_tbl(exp(mstl(log(case.iNTS.ts))))
case.iNTS$case_count_tre <- trend_n$Trend
Seasonal-adjusted typhoid cases
case.typhi.ts <- ts(na.omit(case.typhi$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.typhi.ts+1)))), preserve_index=FALSE)
case.typhi$case_count_sea <-trend_n$value
case.typhi$case_count_sea[case.typhi$case_count_sea < 0] <- 0
setnames(case.typhi, old="case_count", new="case_count_obs")
trend_n <- tk_tbl(exp(mstl(log(case.typhi.ts+1))))
case.typhi$case_count_tre <- trend_n$Trend
Seasonal-adjusted rain
climate.rain.ts <- ts(na.omit(climate.rain$rainfall), frequency=12)
trend_n <- tk_tbl(seasadj(mstl(climate.rain.ts)), preserve_index=FALSE)
climate.rain$rainfall_sea <- trend_n$value
climate.rain$rainfall_sea[climate.rain$rainfall_sea < 0] <- 0
setnames(climate.rain, old="rainfall", new="rainfall_obs")
trend_n <-tk_tbl(mstl(climate.rain.ts))
climate.rain$rain_tre <- trend_n$Trend
Seasonal-adjusted temperature
climate.temp.ts <- ts(na.omit(climate.temp$temperature), frequency=12)
trend_n <-tk_tbl(seasadj(mstl(climate.temp.ts)), preserve_index=FALSE)
climate.temp$temperature_sea <- trend_n$value
climate.temp$temperature_sea[climate.temp$temperature_sea < 0] <- 0
setnames(climate.temp, old="temperature", new="temperature_obs")
trend_n <- tk_tbl(mstl(climate.temp.ts))
climate.temp$temp_tre <- trend_n$Trend

rm(list=ls()[grep("^trend_n", ls())])
Plots of decomposed series
x<-case.iNTS.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="A",xlab="Years (2000-2015)",size=1,colour="orange2",is.date=FALSE)+theme_bw()
y<-case.typhi.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="B",xlab="Years (2000-2015)",size=1,colour="red2",is.date=FALSE)+theme_bw()
z<-climate.rain.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="C",xlab="Years (2000-2015)",size=1,colour="blue2",is.date=FALSE)+theme_bw()
v<-climate.temp.ts %>% mstl() %>% ggfortify:::autoplot.ts(main="D",xlab="Years (2000-2015)",size=1,colour="green2",is.date=FALSE)+theme_bw()
grid.arrange(grobs=list(x,y,z,v), ncol=4, nrow=1)

Plots of seasonal-adjusted series
Th <- theme(plot.title=element_text(hjust=0)) + 
  theme(axis.title.x=element_text(size=10)) + 
  theme(axis.title.y=element_text(size=10)) +
  theme(axis.text.x=element_text(face="bold", size=10), axis.text.y=element_text(face="bold", size=10)) + 
  theme(legend.key.height=unit(1, "line")) + 
  theme(legend.key.width=unit(1, "line")) 

S1 <- ggplot(as.data.frame(case.iNTS)) + 
  geom_line(aes(date, case_count_obs, color="Original"), size=0.8) + 
  geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) + 
  scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="orange2")) +
  labs(title="A", x="", y="NTS cases") + theme(plot.title=element_text(hjust=0)) +
  theme(legend.justification=c(0.5,0), legend.position=c(0.7,0.7), legend.text=element_text(size=10), 
  legend.title=element_text(size=0))

S2<-ggplot(as.data.frame(case.typhi)) + 
  geom_line(aes(date, case_count_obs, color="Original"), size=0.8) + 
  geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) + 
  scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="red2")) +
  labs(title="B", x ="", y="Typhoid cases") + Th +
  theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10), 
  legend.title = element_text(size=0))

S3<-ggplot(as.data.frame(climate.rain)) + 
  geom_line(aes(date, rainfall_obs, color="Original"), size=0.8) + 
  geom_line(aes(date, rainfall_sea, color="Seasonal-adjusted"), size=0.8) + 
  scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="blue2")) +
  labs(title="C", x ="Year", y = "Rainfall (mm)") + Th + 
  theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10), 
  legend.title = element_text(size=0))

S4<-ggplot(as.data.frame(climate.temp)) + 
  geom_line(aes(date, temperature_obs, color="Original"), size=0.8) + 
  geom_line(aes(date, temperature_sea, color="Seasonal-adjusted"), size=0.8) + 
  scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="green2")) + 
  labs(title="D", x ="Year", y = "Temperature (°C)") + Th + 
  theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10), 
  legend.title = element_text(size=0))

grid.arrange(grobs=list(S1, S2, S3, S4), ncol=2, nrow=2)

Inter/extraporate population denominators
census.year <-c(1998, 2008)
census.popn <-c(809397, 1022680)
census.count.1998.2008 <-approx(census.year, census.popn, n=11) 
census.count.2009.2015 <-approxExtrap(census.year, census.popn, xout=c(2009, 2010, 2011, 2012, 2013, 2014, 2015))
Calculate NTS/typhoid incidence
case.iNTS$census[year(case.iNTS$date)==2000]<-852054; case.iNTS$census[year(case.iNTS$date)==2001]<-873382
case.iNTS$census[year(case.iNTS$date)==2002]<-894710; case.iNTS$census[year(case.iNTS$date)==2003]<-916039
case.iNTS$census[year(case.iNTS$date)==2004]<-937367; case.iNTS$census[year(case.iNTS$date)==2005]<-958695
case.iNTS$census[year(case.iNTS$date)==2006]<-980023; case.iNTS$census[year(case.iNTS$date)==2007]<-1001352
case.iNTS$census[year(case.iNTS$date)==2008]<-1022680; case.iNTS$census[year(case.iNTS$date)==2009]<-1044008
case.iNTS$census[year(case.iNTS$date)==2010]<-1065337; case.iNTS$census[year(case.iNTS$date)==2011]<-1086665
case.iNTS$census[year(case.iNTS$date)==2012]<-1107993; case.iNTS$census[year(case.iNTS$date)==2013]<-1129322
case.iNTS$census[year(case.iNTS$date)==2014]<-1150650; case.iNTS$census[year(case.iNTS$date)==2015]<-1171978
case.iNTS$incid_sea <-case.iNTS$case_count_sea*100000/case.iNTS$census 
case.iNTS$incid_obs <-case.iNTS$case_count_obs*100000/case.iNTS$census 

case.typhi$census[year(case.typhi$date)==2000]<-852054; case.typhi$census[year(case.typhi$date)==2001]<-873382
case.typhi$census[year(case.typhi$date)==2002]<-894710; case.typhi$census[year(case.typhi$date)==2003]<-916039
case.typhi$census[year(case.typhi$date)==2004]<-937367; case.typhi$census[year(case.typhi$date)==2005]<-958695
case.typhi$census[year(case.typhi$date)==2006]<-980023; case.typhi$census[year(case.typhi$date)==2007]<-1001352
case.typhi$census[year(case.typhi$date)==2008]<-1022680; case.typhi$census[year(case.typhi$date)==2009]<-1044008
case.typhi$census[year(case.typhi$date)==2010]<-1065337; case.typhi$census[year(case.typhi$date)==2011]<-1086665
case.typhi$census[year(case.typhi$date)==2012]<-1107993; case.typhi$census[year(case.typhi$date)==2013]<-1129322
case.typhi$census[year(case.typhi$date)==2014]<-1150650; case.typhi$census[year(case.typhi$date)==2015]<-1171978
case.typhi$incid_sea <-case.typhi$case_count_sea*100000/case.typhi$census 
case.typhi$incid_obs <-case.typhi$case_count_obs*100000/case.typhi$census
Plots of Age-sex distribution of cases
case$sex[case$sex == ""] <- NA
case$age[case$age == ""] <- NA
case$date <- ymd(case$case_date)
case$year <- year(case$date)

agesex.p1 <- ggplot(subset(case, organism=="iNTS" & !is.na(age) & sex != "Unknown"), aes(x=age, color=sex)) + 
  geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) + 
  scale_color_manual(values=c(Male="gray40",Female="red3")) + 
  theme_bw() + 
  scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) + 
  labs(title="A", x="Age (years)", y="Number of NTS cases") + 
  theme(plot.title=element_text(hjust=0,face="bold")) +
  theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) + 
  theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) + 
  labs(color="Missing sex: 2,596 (32.3%)") + 
  theme(legend.key.height=unit(1,"line")) + 
  theme(legend.key.width=unit(1,"line"))

agesex.p2<- ggplot(subset(case, organism=="typhi" & !is.na(age) & sex !="Unknown"), aes(x=age, color=sex)) + 
  geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) + 
  scale_color_manual(values=c(Male="gray40",Female="red3")) + 
  theme_bw() + 
  scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) + 
  labs(title="B", x="Age (years)", y="Number of typhoid cases") + 
  theme(plot.title = element_text(hjust=0,face="bold")) +
  theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) + 
  theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) + 
  labs(color="Missing sex: 61 (2.4%)") + 
  theme(legend.key.height=unit(1,"line")) + 
  theme(legend.key.width=unit(1,"line"))

grid.arrange(grobs=list(agesex.p1, agesex.p2), ncol=2, nrow=1)

Generating yearly and monthly dynamics
case.iNTS.spi <- subset(case.iNTS, year(case.iNTS$date)<2011, select=c(date,incid_sea)) 
case.iNTS.spi$month <-month(case.iNTS.spi$date)
case.iNTS.spi$year <- year(case.iNTS.spi$date); case.iNTS.spi$date <- NULL
case.iNTS.spi <- spread(case.iNTS.spi, year, incid_sea); case.iNTS.spi <- as.matrix(case.iNTS.spi)[,-1]
YMD1 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~case.iNTS.spi, type="contour", colorscale='jet', contours=list(showlabels=TRUE)) %>% 
colorbar(title="NTS incidence per \n 100,000 population") %>%
layout(title="A", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
case.typhi.spi <- subset(case.typhi, year(case.typhi$date)>2010, select=c(date,incid_sea)) 
case.typhi.spi$month <- month(case.typhi.spi$date)
case.typhi.spi$year <- year(case.typhi.spi$date); case.typhi.spi$date <- NULL
case.typhi.spi <- spread(case.typhi.spi, year, incid_sea); case.typhi.spi <- as.matrix(case.typhi.spi)[,-1]
YMD2 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~case.typhi.spi, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% 
colorbar(title="Typhoid incidence per \n 100,000 population") %>%
layout(title="D", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spin <- subset(climate.rain, year(climate.rain$date)<2011, select=c(date,rainfall_obs)) 
climate.rain.spin$month <- month(climate.rain.spin$date)
climate.rain.spin$year <- year(climate.rain.spin$date); climate.rain.spin$date <- NULL
climate.rain.spin <- spread(climate.rain.spin, year, rainfall_obs); climate.rain.spin <- as.matrix(climate.rain.spin)[,-1]
YMD3 <- plot_ly(x = c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.rain.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% colorbar(title="Rainfall (mm)") %>% 
layout(title="B", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spit <- subset(climate.rain, year(climate.rain$date)>2010, select=c(date,rainfall_obs)) 
climate.rain.spit$month <- month(climate.rain.spit$date)
climate.rain.spit$year <- year(climate.rain.spit$date); climate.rain.spit$date <- NULL
climate.rain.spit <- spread(climate.rain.spit, year, rainfall_obs); climate.rain.spit <- as.matrix(climate.rain.spit)[,-1]
YMD4 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.rain.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% 
colorbar(title="Rainfall (mm)") %>%
layout(title="E", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spin <- subset(climate.temp, year(climate.temp$date)<2011, select=c(date,temperature_obs)) 
climate.temp.spin$month <- month(climate.temp.spin$date)
climate.temp.spin$year <- year(climate.temp.spin$date); climate.temp.spin$date <- NULL
climate.temp.spin <- spread(climate.temp.spin, year, temperature_obs); climate.temp.spin <- as.matrix(climate.temp.spin)[,-1]
YMD5 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.temp.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% 
colorbar(title="Temperature (°C)") %>%
layout(title="C", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spit <- subset(climate.temp, year(climate.temp$date)>2010, select=c(date,temperature_obs)) 
climate.temp.spit$month <- month(climate.temp.spit$date)
climate.temp.spit$year <- year(climate.temp.spit$date); climate.temp.spit$date <- NULL
climate.temp.spit <- spread(climate.temp.spit, year, temperature_obs); climate.temp.spit <- as.matrix(climate.temp.spit)[,-1]
YMD6 <- plot_ly(x=c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"), 
z=~climate.temp.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Temperature (°C)") %>%
layout(title="F", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size = 13))
Plots of yearly and monthly dynamics
htmltools::tags$div(
  style = "display: flex; flex-wrap: wrap",
  tags$div(YMD1, style = "width: 33%; padding: 1em;"),
  tags$div(YMD3, style = "width: 33%; padding: 1em;"),
  tags$div(YMD5, style = "width: 33%; padding: 1em;"),
  tags$div(YMD2, style = "width: 33%; padding: 1em;"),
  tags$div(YMD4, style = "width: 33%; padding: 1em;"),
  tags$div(YMD6, style = "width: 33%; padding: 1em;")
)
Constructing DLNM datasets for NTS and Typhoid
mo.dlnmN <- bind_cols(case.iNTS, climate.rain, climate.temp, id=NULL)
mo.dlnmN$date1 <- mo.dlnmN$date2 <- NULL
mo.dlnmN <- subset(mo.dlnmN, year(date) < 2011)
mo.dlnmN$time <- seq.int(from = 1, to=132, by=1)
mo.dlnmN$year <- year(mo.dlnmN$date)
mo.dlnmN$month <- month(mo.dlnmN$date)
mo.dlnmN$incid_seaX<-round(mo.dlnmN$incid_sea, digits = 0)
mo.dlnmN$incid_obsX<-round(mo.dlnmN$incid_obs, digits = 0)

mo.dlnmT <- bind_cols(case.typhi, climate.rain, climate.temp, id=NULL)
mo.dlnmT$date1 <- mo.dlnmT$date2 <- NULL
mo.dlnmT <- subset(mo.dlnmT, year(date) > 2010)
mo.dlnmT$time <- seq.int(from = 1, to=60, by=1)
mo.dlnmT$year <- year(mo.dlnmT$date)
mo.dlnmT$month <- month(mo.dlnmT$date)
mo.dlnmT$incid_seaX<-round(mo.dlnmT$incid_sea, digits = 0)
AIC to QAIC for NTS and Typhoid
nts.quasipoisson <- function(...) { 
  res <- quasipoisson(...)
  res$aic <- poisson(...)$aic 
  res
}
typ.quasipoisson <- function(...) { 
  res <- quasipoisson(...) 
  res$aic <- poisson(...)$aic 
  res
}
Calculate QAIC values to select optimal models
QAICtable <- data.frame(model.no=rep(NA,27), lag.df=rep(NA,27), fx1.df=rep(NA,27), fx2.df=rep(NA,27), 
                      QAIC.ntsR=rep(NA,27), QAIC.ntsT=rep(NA,27), QAIC.nts=rep(NA,27), QAIC.typR=rep(NA,27), 
                      QAIC.typT=rep(NA,27), QAIC.typ=rep(NA,27))
l=1
for(k in 3:5){
  for(j in 3:5){
    for(i in 3:5){
  nts.lagknots <- logknots(8, fun="ns", df=i)
  nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=j)
  nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=k)
  nts.mo.cb.rain <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots))
  nts.mo.cb.temp <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots))
  nts.modelR <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
  nts.modelT <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.temp + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
  nts.model <-  glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + nts.mo.cb.temp + year, family = nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
  
  typ.lagknots <- logknots(8, fun="ns", df=i)
  typ.varknots.r=equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=j)
  typ.varknots.t=equalknots(mo.dlnmT$temperature_obs, fun="ns", df=k)
  typ.mo.cb.rain <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r), arglag=list(knots=typ.lagknots))
  typ.mo.cb.temp <- crossbasis(mo.dlnmT$temperature_obs, lag =8, argvar = list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
  typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
  typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
  typ.model <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
  
  QAICtable[l,] <- c(l,i,j,k, QAIC(nts.modelR, chat=summary(nts.modelR)$dispersion), QAIC(nts.modelT, chat=summary(nts.modelT)$dispersion), 
                     QAIC(nts.model, chat=summary(nts.model)$dispersion), QAIC(typ.modelR, chat=summary(typ.modelR)$dispersion), QAIC(typ.modelT, chat=summary(typ.modelT)$dispersion), 
                     QAIC(typ.model, chat=1))
  l=l+1  
  }
  }
}
kable(QAICtable)
model.no lag.df fx1.df fx2.df QAIC.ntsR QAIC.ntsT QAIC.nts QAIC.typR QAIC.typT QAIC.typ
1 3 3 3 570.8194 565.8631 562.3950 218.5004 220.2006 212.6658
2 4 3 3 576.3831 571.4838 573.3630 223.1714 216.1853 220.6390
3 5 3 3 581.4852 574.9872 583.4617 228.2399 221.1698 230.8468
4 3 4 3 571.1280 565.8631 565.3731 211.3436 220.2006 216.2219
5 4 4 3 578.3398 571.4838 578.1025 218.1843 216.1853 228.0142
6 5 4 3 585.4681 574.9872 590.3358 222.7394 221.1698 237.6362
7 3 5 3 572.7386 565.8631 566.1964 216.3176 220.2006 221.5570
8 4 5 3 581.6967 571.4838 581.3661 222.4327 216.1853 234.0871
9 5 5 3 590.5606 574.9872 594.8792 227.8628 221.1698 246.2023
10 3 3 4 570.8194 571.7243 566.2250 218.5004 220.7893 215.5219
11 4 3 4 576.3831 577.5401 578.9833 223.1714 220.2852 225.9024
12 5 3 4 581.4852 583.4447 589.9669 228.2399 226.2618 236.0380
13 3 4 4 571.1280 571.7243 568.6740 211.3436 220.7893 220.0702
14 4 4 4 578.3398 577.5401 583.3435 218.1843 220.2852 233.1946
15 5 4 4 585.4681 583.4447 597.5272 222.7394 226.2618 245.6773
16 3 5 4 572.7386 571.7243 571.2271 216.3176 220.7893 224.5520
17 4 5 4 581.6967 577.5401 588.2408 222.4327 220.2852 239.9495
18 5 5 4 590.5606 583.4447 603.1850 227.8628 226.2618 252.6551
19 3 3 5 570.8194 574.0369 566.7772 218.5004 220.7423 220.4298
20 4 3 5 576.3831 580.7484 580.7418 223.1714 220.9285 232.7209
21 5 3 5 581.4852 587.7540 593.8723 228.2399 225.5709 245.1746
22 3 4 5 571.1280 574.0369 569.9469 211.3436 220.7423 224.8192
23 4 4 5 578.3398 580.7484 586.2703 218.1843 220.9285 239.8864
24 5 4 5 585.4681 587.7540 601.6447 222.7394 225.5709 254.6679
25 3 5 5 572.7386 574.0369 575.0626 216.3176 220.7423 229.8974
26 4 5 5 581.6967 580.7484 593.7784 222.4327 220.9285 246.7405
27 5 5 5 590.5606 587.7540 610.6531 227.8628 225.5709 261.6885
Reformulate cross-basis matrix using selected best NTS model
sort(mo.dlnmN$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmN$temperature_obs, decreasing=FALSE)
nts.lagknots.r <- logknots(8, fun="ns", df=3)
nts.varknots.r <- equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=3)
nts.lagknots.t <- logknots(8, fun="ns", df=3)
nts.varknots.t <- equalknots(mo.dlnmN$temperature_obs, fun="ns", df=3)
nts.cb.rain <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots.r))
nts.cb.temp <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots.t))
summary(nts.cb.rain)
summary(nts.cb.temp)
Fit NTS model
nts.model <- glm(mo.dlnmN$incid_seaX ~ nts.cb.rain + nts.cb.temp + year, family=quasipoisson(), na.action=na.delete, mo.dlnmN)
Investigate deviance residuals for a fitted NTS model
par(mfrow=c(1,2))
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Partial autocorrelation (original modelR+T)", xlim=c(0,8))
nts.model <- update(nts.model,.~.+Lag(residuals(nts.model,type="deviance"), 1)) #add residuals at lag 1 
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Partial autocorrelation (adjusted modelR+T)", xlim=c(0,8))

Run predictions
nts.pred.rain <- crosspred(nts.cb.rain, nts.model, cen=0, by=0.2)
nts.pred.temp <- crosspred(nts.cb.temp, nts.model, cen=23, by=0.2)
Plots of NTS predictions due to rainfall
plot(nts.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1.1, 
    cex.axis=1.1, col="gray80",main="A")

plot(nts.pred.rain, "contour", key.title=title("NTS"), plot.title=title("B", xlab ="Rainfall (mm)", ylab = "Lag (month)", 
    cex.lab=1.1, cex.axis=1.1))

plot(nts.pred.rain, "slices", xlab="Rainfall (mm)", lag=c(2), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)), 
    ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="C")

plot(nts.pred.rain, "slices", xlab="Lag (month)", var=c(9), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)), 
    ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="D")

Plots of NTS predictions due to temperature
plot(nts.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1, cex.axis=1, col="gray80",main="A")

plot(nts.pred.temp, "contour", key.title=title("NTS"), plot.title=title("B", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1, cex.axis=1))

plot(nts.pred.temp, xlab="Temperature (°C)", "slices",lag=c(3), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="C")

plot(nts.pred.temp, xlab="Lag (month)", "slices",var=c(19), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.1, cex.axis=1.1,main="D")

Reformulate cross-basis matrix using selected best Typhoid models
sort(mo.dlnmT$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmT$temperature_obs, decreasing=FALSE)
typ.lagknots <- logknots(8, fun="ns", df=3)
typ.varknots.r1 <- equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=4)
typ.varknots.r2 <- equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=3)
typ.varknots.t <- equalknots(mo.dlnmT$temperature_obs, fun="ns", df=3)
typ.cb.rain1 <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r1), arglag=list(knots=typ.lagknots))
typ.cb.rain2 <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r2), arglag=list(knots=typ.lagknots))
typ.cb.temp <- crossbasis(mo.dlnmT$temperature_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
summary(typ.cb.rain1)
summary(typ.cb.rain2)
summary(typ.cb.temp)
Fit Typhoid model
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.cb.rain1 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.cb.temp + typ.cb.rain2 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
Investigate deviance residuals for fitted tyhpoid models
par(mfrow=c(1,3))
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="Partial Autocorrelation (original modelR)",xlim=c(0,8))
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Partial Autocorrelation (original modelT)",xlim=c(0,8))
typ.modelT <- update(typ.modelT,.~.+Lag(residuals(typ.modelT,type="deviance"),1)) #add residuals at lag 1
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Partial Autocorrelation (adjusted modelT)",xlim=c(0,8))

Run predictions
typ.pred.rain <- crosspred(typ.cb.rain1, typ.modelR, cen=0, by=0.2)
typ.pred.temp <- crosspred(typ.cb.temp, typ.modelT, cen=23, by=0.2)
Plots of TYP predictions due to rainfall
plot(typ.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of typhoid", theta=150, phi=5, lphi=150, col="gray80", main="A", cex.lab=1.1, cex.axis=1.1)

plot(typ.pred.rain, "contour", key.title=title("TYP"), plot.title=title("B", xlab ="Rainfall (mm)", ylab = "Lag (month)", cex.lab=1.1, cex.axis=1.1))

plot(typ.pred.rain, xlab="Rainfall (mm)", "slices", lag=c(4), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="C")

plot(typ.pred.rain, xlab="Lag (month)", "slices", var=c(9), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="D")

Plots of TYP predictions due to temperature
plot(typ.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of typhoid", theta=150, phi=5, lphi=100, col="gray80", main="A", cex.lab=1.1, cex.axis=1.1)

plot(typ.pred.temp, "contour", key.title=title("TYP"), plot.title=title("B", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1.1, cex.axis=1.1))

plot(typ.pred.temp, xlab="Temperature (°C)", "slices", lag=c(2), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="C")

plot(typ.pred.temp, xlab="Temperature (°C)", "slices", lag=c(5), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="D")

plot(typ.pred.temp, xlab="Lag (month)", "slices", var=c(19), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="E")

plot(typ.pred.temp, xlab="Lag (month)", "slices", var=c(25), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.1, cex.axis=1.1,main="F")

Plots of validation checks
par(mfrow=c(3,4))
plot(mo.dlnmN$date, residuals(nts.model,type="deviance"), pch=19, cex=0.8, col=grey(0.4),main="A", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2) 
abline(h=0,lty=2,lwd=3)
pacf(residuals(nts.model,type="deviance"),na.action=na.omit,main="B",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
acf(residuals(nts.model,type="deviance"),na.action=na.omit,main="C",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
plot(mo.dlnmN$incid_seaX, pch=10, cex=0.8, col=grey(0.6), size =1 ,main="D", ylab="NTS incidence",xlab="",cex.lab=1.2, cex.axis=1.2)
lines(predict(nts.model, type="response"), col="orange2",lwd=3)
legend("topright", legend=c("observed", "predicted"), col=c("grey", "orange2"), lty=1:1, cex=0.8, lwd=3)

plot(mo.dlnmT$date, residuals(typ.modelR,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="E", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2) 
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="F",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="G",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
plot(mo.dlnmT$incid_seaX, pch=10, cex=0.8, col=grey(0.6), size =1 ,main="H", ylab="Typhoid incidence",xlab="",cex.lab=1.2, cex.axis=1.2)
lines(predict(typ.modelR, type="response"), col="red",lwd=3)
legend("topleft", legend=c("observed", "predicted"), col=c("grey", "red2"), lty=1:1, cex=0.8, lwd=3)

plot(mo.dlnmT$date, residuals(typ.modelT,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="I", ylab="Residuals", xlab="Year",cex.lab=1.2,cex.axis=1.2) 
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="J",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="K",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
plot(mo.dlnmT$incid_seaX, pch=10, cex=0.8, col=grey(0.6), size =1 ,main="L", ylab="Typhoid incidence", xlab="Month",cex.lab=1.2, cex.axis=1.2)
lines(predict(typ.modelT, type="response"), col="red2",lwd=3)
legend("topleft", legend=c("observed", "predicted"), col=c("grey", "red"), lty=1:1, cex=0.8, lwd=3)

Plots of TYP predictions due to temperature
cbind(nts.pred.rain$matRRfit, nts.pred.rain$matRRlow, nts.pred.rain$matRRhigh)["9",]
##      lag0      lag1      lag2      lag3      lag4      lag5      lag6 
## 1.2248467 1.2514626 1.2554780 1.2238951 1.1635420 1.0841819 0.9951394 
##      lag7      lag8      lag0      lag1      lag2      lag3      lag4 
## 0.9042881 0.8176177 1.0254867 1.1162267 1.1453653 1.1143857 1.0596817 
##      lag5      lag6      lag7      lag8      lag0      lag1      lag2 
## 0.9906237 0.9027750 0.7975713 0.6897005 1.4629635 1.4030831 1.3761766 
##      lag3      lag4      lag5      lag6      lag7      lag8 
## 1.3441658 1.2775817 1.1865762 1.0969538 1.0252839 0.9692593
cbind(nts.pred.temp$matRRfit, nts.pred.temp$matRRlow, nts.pred.temp$matRRhigh)["19",]
##      lag0      lag1      lag2      lag3      lag4      lag5      lag6 
## 1.0822001 1.1309828 1.1629884 1.1659357 1.1432279 1.1012306 1.0467408 
##      lag7      lag8      lag0      lag1      lag2      lag3      lag4 
## 0.9861515 0.9249535 0.9602616 1.0318836 1.0674043 1.0690932 1.0514755 
##      lag5      lag6      lag7      lag8      lag0      lag1      lag2 
## 1.0192753 0.9708904 0.9049599 0.8287981 1.2196229 1.2395991 1.2671318 
##      lag3      lag4      lag5      lag6      lag7      lag8 
## 1.2715505 1.2429867 1.1897756 1.1285169 1.0746274 1.0322647
cbind(nts.pred.temp$matRRfit, nts.pred.temp$matRRlow, nts.pred.temp$matRRhigh)["25",]
##      lag0      lag1      lag2      lag3      lag4      lag5      lag6 
## 0.9983030 0.9778710 0.9641177 0.9603166 0.9651159 0.9768977 0.9941390 
##      lag7      lag8      lag0      lag1      lag2      lag3      lag4 
## 1.0153070 1.0387806 0.9057049 0.9118456 0.9084918 0.9062719 0.9129857 
##      lag5      lag6      lag7      lag8      lag0      lag1      lag2 
## 0.9274517 0.9428357 0.9521165 0.9547300 1.1003682 1.0486771 1.0231495 
##      lag3      lag4      lag5      lag6      lag7      lag8 
## 1.0175842 1.0202227 1.0289798 1.0482338 1.0826914 1.1302307
cbind(typ.pred.rain$matRRfit, typ.pred.rain$matRRlow, typ.pred.rain$matRRhigh)["9",]
##      lag0      lag1      lag2      lag3      lag4      lag5      lag6 
## 0.9870340 1.2979004 1.6063110 1.8078055 1.8722857 1.8143089 1.6725894 
##      lag7      lag8      lag0      lag1      lag2      lag3      lag4 
## 1.4915122 1.3081090 0.6284634 0.9438326 1.2552773 1.4434731 1.5069006 
##      lag5      lag6      lag7      lag8      lag0      lag1      lag2 
## 1.4634299 1.3249680 1.1234975 0.9124989 1.5501876 1.7847926 2.0555100 
##      lag3      lag4      lag5      lag6      lag7      lag8 
## 2.2640954 2.3262675 2.2493164 2.1114137 1.9800744 1.8752341
cbind(typ.pred.temp$matRRfit, typ.pred.temp$matRRlow, typ.pred.temp$matRRhigh)["19",]
##      lag0      lag1      lag2      lag3      lag4      lag5      lag6 
## 1.4960476 1.5619360 1.5874626 1.5468052 1.4526163 1.3245023 1.1812571 
##      lag7      lag8      lag0      lag1      lag2      lag3      lag4 
## 1.0380766 0.9055470 1.0424914 1.1457962 1.0968536 1.0279449 0.9761052 
##      lag5      lag6      lag7      lag8      lag0      lag1      lag2 
## 0.9252186 0.8373032 0.6896236 0.5229565 2.1469324 2.1292129 2.2975149 
##      lag3      lag4      lag5      lag6      lag7      lag8 
## 2.3275628 2.1617487 1.8960993 1.6665030 1.5625959 1.5680377
cbind(typ.pred.temp$matRRfit, typ.pred.temp$matRRlow, typ.pred.temp$matRRhigh)["25",]
##      lag0      lag1      lag2      lag3      lag4      lag5      lag6 
## 0.6881750 0.9207072 1.1640344 1.3467161 1.4416997 1.4504554 1.3928588 
##      lag7      lag8      lag0      lag1      lag2      lag3      lag4 
## 1.2966573 1.1885051 0.4844865 0.7123388 0.8891649 0.9975957 1.0683040 
##      lag5      lag6      lag7      lag8      lag0      lag1      lag2 
## 1.1081654 1.1070709 1.0400232 0.9069586 0.9774984 1.1900260 1.5238750 
##      lag3      lag4      lag5      lag6      lag7      lag8 
## 1.8180153 1.9456054 1.8984720 1.7524222 1.6166181 1.5574519